Method for solving high-dimensional nonlinear filtering problem

ABSTRACT

A method for solving high-dimensional nonlinear filtering problems is revealed. The method uses a fast computational module to solve an equation and get approximate numerical solutions of a signal-observation model. The equation-solving process of the fast computational module is speeded up by a transformation module and the computational stability is further improved. D-dimensional nonlinear filtering problems are solved and approximate numerical solutions are obtained based on Yau-Yau filtering theory. A Quasi-Implicit Euler Method (QIEM) is applied to solve the Kolmogorov equations and estimate approximate numerical solutions of the signal-observation model. Moreover, QIEM is more efficient and the numerical solutions are more stable by Fast Fourier transformation (FFT) acceleration.

BACKGROUND OF THE INVENTION

1. Field of the invention

The present invention relates to a method for solving high-dimensional nonlinear filtering problems, especially to a method for solving high-dimensional nonlinear filtering problems that solves D-dimensional nonlinear filtering problems and gets approximate numerical solutions of a signal-observation model based on Yau-Yau filtering theory. A Quasi-Implicit Euler Method (QIEM) is applied to solve the Kolmogorov equations and estimate approximate states of a given signal-observation model. By Fast Fourier transformation (FFT) acceleration, QIEM is more efficient and the numerical solutions are more stable.

2. Description of Related Art

The nonlinear filtering problem has a variety of applications in military, engineering and commercial industries. The core issue of the nonlinear filtering problem is to solve the Duncan-Mortensen-Zakai (DMZ) equation in real time. Yau and Yau have proved that the real-time solution of the DMZ equation can be reduced to an off-time solution of the Kolmogorov equation. Based on the work of Yau and Yau, Liu et al. proposed as numerical method to solve the nonlinear filtering problem by explicit finite difference schemes (Existence and uniqueness and decay estimates for the time dependent parabolic equation with application to Duncan-Mortensen-Zakai equation, Asian Journal of Mathematics, 2:1079-1149, 1998). In order to improve the stability of the algorithm proposed by Liu, an efficient and reliable quasi-implicit numerical scheme is proposed for solving the Kolmogorov equations and estimate approximate states of a given signal-observation model.

SUMMARY OF THE INVENTION

Therefore it is a primary object of the present invention to provide a method for solving high-dimensional nonlinear filtering problems by which high-dimensional nonlinear filtering problems are solve and approximate numerical solutions of a signal-observation model are obtained based on Yau-Yau filtering theory. A Quasi-Implicit Euler Method (QIEM) is applied to solve Kolmogorov equations and estimate approximate numerical solutions of the signal-observation model. Moreover, stability of the numerical solutions is ensured by fast Fourier transformation (FFT) applied in the QIEM.

In order to achieve the above object, a method for solving high-dimensional nonlinear filtering problems of the present invention solves an equation and obtains approximate numerical solutions of a signal-observation model by using a fast computational module. Moreover, the equation-solving process of the fast computational module is speeded up by a transformation module in order to improve the computational stability.

In the fast computational module, a Quasi-Implicit Euler Method (QIEM) is applied to solve Kolmogorov equations and estimate approximate numerical solutions of the given signal-observation model.

In the method mentioned above, the Quasi-Implicit Euler Method(QIEM) is iteratively formulated by:

${{\left\lbrack {I_{N^{D}} - {\frac{\Delta t}{2\left( {\Delta \; s} \right)^{2}}L_{N}^{(D)}}} \right\rbrack u_{n + 1}} = {\left\lbrack {I_{N^{D}} + {{\Delta t}\left( {{\frac{1}{2\Delta \; s}K_{N}^{(D)}} + Q_{N}^{(D)}} \right)}} \right\rbrack u_{N}}};$ ${{{wherein}\mspace{14mu} L_{N}^{(D)}} = {\sum\limits_{d = 1}^{D}\; \left\lbrack {I_{N^{D - d}} \otimes L_{N} \otimes I_{N^{d - 1}}} \right\rbrack}},{{K_{N}^{(D)} = {\sum\limits_{d = 1}^{D}{P_{d}\;\left\lbrack {I_{N^{D - d}} \otimes K_{N} \otimes I_{N^{d - 1}}} \right\rbrack}}};}$

-   -   wherein I_(N) is the identity matrix of size N and         P_(d)=diag{p_(d)(s_(j))}_(j=1) ^(N) ^(D) ,         Q=diag{q(s_(j))}_(j=1) ^(N) ^(D) ,     -   wherein the matrices L_(N) and K_(N) are defined by

${L_{N} = \begin{bmatrix} {- 2} & 1 & \; & \; \\ 1 & {- 2} & \ldots & \; \\ \; & \ldots & \ldots & 1 \\ \; & \; & 1 & {- 2} \end{bmatrix}},{K_{N} = \begin{bmatrix} 0 & 1 & \; & \; \\ {- 1} & 0 & \ldots & \; \\ \; & \ldots & \ldots & 1 \\ \; & \; & {- 1} & 0 \end{bmatrix}}$

Thereby the method of the present invention solves D-dimensional nonlinear filtering problems and gets approximate numerical solutions of the signal-observation model based on Yau-Yau filtering theory. The approximate numerical solutions of the signal-observation model are obtained by applying Quasi-Implicit Euler Method (QIEM) to solve the Kolmogorov equations. The QIEM is more efficient by fast Fourier transformation (FFT) acceleration. Thus stability of the numerical solutions is ensured and computational cost is significantly saved. Moreover, the numerical solution of the Kolmogorov equations in each iteration is nonnegative. The probability density functions are preserved in the iterative process. The results prove that the present method has high efficiency and great potential.

DETAILED DESCRIPTION OF THE PREFFERED EMBODIMENT

In order to learn features and functions of the present invention, please refer to the following embodiments with details.

A method for solving high-dimensional nonlinear filtering problems of the present invention solves an equation and gets approximate numerical solutions of a signal-observation model by using a fast computational module. A transformation module is used to accelerate the equation-solving process of the fast computational module for improving the computational stability. In the fast computational module, a Quasi-Implicit Euler Method (QIEM) is applied to solve the Kolmogorov equations and estimate approximate numerical solutions of the signal-observation model.

The nonlinear filtering problem considered here is to determine approximate states for a given observation history of the following signal-observation model:

$\begin{matrix} \left\{ \begin{matrix} {{{dX}(t)} = {{{f\left( {X(t)} \right)}{dt}} + {{dv}(t)}}} \\ {{{dY}(t)} = {{{h\left( {X(t)} \right)}{dt}} + {{dw}(t)}}} \end{matrix} \right. & (1) \end{matrix}$

wherein X(0)=X₀, Y(0)=0, and X(t)=(x_(l)(t), . . . , x_(D)(t))^(T) ∈ R^(D), Y(t)=(y₁(t), . . . , y_(M)(t))^(T) ∈ R^(M) are the state and the measurement/observation vectors at time t, respectively. ƒ(X)=(ƒ₁(X), . . . , ƒ_(D)(X))^(T), h(X)=(h₁(X), . . . , h_(M)(X))^(T) are given vector-valued functions, ν ∈ R^(D) and w ∈ R^(M) are mutually independent standard Brownian processes. From the main results of Yau and Yau, the sate vector x(t) can be estimated from the observation vectors {y(s)|s ∈ [0,t]} by solving the Kolmogorov equations. Specifically, suppose that a set of observations {y(τ₀), . . . , y(τ_(N), )} is measured. For each time period [τ_(k−1), τ_(k)], k=1, . . . , N_(t), the Kolmogorov equations of the from are solved,

$\begin{matrix} \left\{ \begin{matrix} {{{\frac{\partial{\overset{\sim}{u}}_{k}}{\partial t}\left( {t,s} \right)} = {{\frac{1}{2}\Delta {{\overset{\sim}{u}}_{k}\left( {t,s} \right)}} + {\sum\limits_{d = 1}^{D}\; {{{Pd}(s)}\frac{\partial{\overset{\sim}{u}}_{k}}{\partial s_{d}}\left( {t,s} \right)}} + {{q(s)}{{\overset{\sim}{u}}_{k}\left( {t,s} \right)}}}},{t \in \left\lbrack {\tau_{k - 1},\tau_{k}} \right\rbrack},} \\ {{{{\overset{\sim}{u}}_{k}\left( {\tau_{k - 1},s} \right)} = {\exp \left\{ {\sum\limits_{j = 1}^{M}\; {\left\lbrack {{y_{j}\left( \tau_{k - 1} \right)} - {y_{j}\left( \tau_{k - 2} \right)}} \right\rbrack {h_{j}(s)}}} \right\} {{\overset{\sim}{u}}_{k - 1}\left( {\tau_{k - 1},s} \right)}}},} \\ {{{{\overset{\sim}{u}}_{1}\left( {0,s} \right)} = {{\sigma_{0}(s)}\exp \left\{ {\sum\limits_{j = 1}^{M}\; {{y_{j}\left( \tau_{0} \right)}{h_{j}(s)}}} \right\}}},} \end{matrix} \right. & (2) \\ \begin{matrix} {\mspace{79mu} {{{{for}\mspace{14mu} k} = 2},\ldots \mspace{14mu},{{N_{\tau}\mspace{14mu} {and}\mspace{14mu} \tau} = 0},{{{wherein}\mspace{14mu} \Delta} = {\sum\limits_{i = 1}^{D}\; \frac{\partial^{2}}{\partial s_{i}^{2}}}},}} & \; \end{matrix} & \; \\ {\mspace{79mu} {{{p_{d}(s)} = {- {f_{d}(s)}}},{d = 1},\ldots \mspace{14mu},D,}} & (3) \\ {\mspace{79mu} {{{{and}\mspace{14mu} q(s)} = {- \left\lbrack {{\sum\limits_{d = 1}^{D}\; {\frac{\partial f_{d}}{\partial s_{d}}(s)}} + {\frac{1}{2}{\sum\limits_{j = 1}^{M}\; {h_{j}^{2}(s)}}}} \right\rbrack}},}} & (4) \end{matrix}$

for t ∈[τ_(k−1), τ_(k)], k=2, . . . , N, the expectation of ũ_(k)(t,s) is computed with respect to s_(d)over R^(D) by

$\begin{matrix} {{{{\hat{x}}_{d}(t)} = {\int_{R^{D}}^{\;}{s_{d}{\overset{\sim}{u}\ }_{k}\left( {t,s} \right){s}}}},} & (5) \end{matrix}$

for d=1, . . . , D. Then the state vector X(t) can be estimated by {circumflex over (X)}(t)=({circumflex over (x)}₁(t), . . . , {circumflex over (x)}_(D)(t))^(T).

In order to solve the nonlinear filtering problem (1), the states and observations {X_(k), Y_(k)}_(k=0) ^(N) ^(r) by using Euler forward difference method with Gaussian noise are first generated. Based on the Yau-Yau method, an implicit Euler method (IEM) for solving the Kolmogorov equation (2) which is stable and reliable, but costly is first proposed. Furthermore, the quasi-implicit Euler method (QIEM) is developed for solving the Kolmogorov equation (2) which is also stable and reliable, but much more efficient because the fast Fourier transformation (FFT) can be applied in the QIEM.

Given a terminal time Γ, a set of states and observations {X_(k), Y_(k)}_(k=0) ^(N) ^(r) are generated by Euler forward difference method. A time interval [0, Γ] is partitioned uniformly as

P _([0, Γ])={0=τ₀<τ₁< . . . <τ_(N), =Γ},   (6)

wherein τ_(k)−τ_(k−1)=Δτ, k=1, 2, . . . , N_(r). With the Euler forward discretization, the signal observation model in (1) an be stimulated by

$\quad\left\{ \begin{matrix} {{X_{k + 1} = {X_{k} + {{f\left( X_{k} \right)}{\Delta\tau}} + {v\sqrt{\Delta\tau}}}},} \\ {{Y_{k + 1} = {Y_{k} + {{h\left( X_{k} \right)}{\Delta\tau}} + {w\sqrt{\Delta\tau}}}},} \end{matrix} \right.$

wherein X₀ is the initial vector and Y₀ is zero, Δτ is the size of the time step, ν and w are mutually independent Brownian motion with ν,w˜N(0,1). The algorithm of “States and Observation Generator” is stated as the following Algorithm 1.

Input: a terminal time Γ, time step Δτ, the initial state X₀, the vector-valued functions f, h Output: the state X(t) and observation Y(t) at t = τ₀, . . . , τ_(N),  1. $N_{\tau} = {\frac{\Gamma}{\Delta\tau} + 1.}$  2. Set the initial state X(τ₀) = X₀.  3. Generate v,w ~ N(0,1), v ⊥ w.  4. for k = 1 to N_(τ) do  5. X(τ_(k)) = X(τ_(k−1)) + f(X(τ_(k−1)))Δτ + v{square root over (Δτ)}.  6. end for  7. Y(τ₀) = 0.  8. for k = 1 to N_(τ) do  9. Y(τ_(k)) = X(τ_(k−1)) + h(X(τ_(k−1)))Δτ + w{square root over (Δτ.)} 10. end for

The IEM is proposed for solving the Kolmogorov equations (2). From the equation (5), the state vector X(t) can be estimated by the solution of the equation (2). For each time interval, [τ_(k−1), τ_(k)], k=1, . . . , N_(r), is partitioned uniformly by

P _([r) _(k−1) _(, r) _(t) _(])={τ_(k−1) =t ₀ ^((k)) <t ₁ ^((k)) <. . . <t _(N),^((k))=τ_(k)},

wherein t_(n) ^((k)−t) _(n−1) ^((k)=Δt, n=)1, . . . , N_(r). Then the partition

$P_{\lbrack{0,\Gamma}\rbrack}^{*} = {{\overset{N_{\tau}}{\bigcup\limits_{k = 1}}P_{\lbrack{\tau_{k - 1},\tau_{k}}\rbrack}} = \left\{ {0 = {\tau_{0} = {{t_{0}^{(1)} < \ldots < t_{N_{t}}^{(1)}} = {\tau_{1} = {{t_{0}^{(2)} < \ldots < t_{N_{t}}^{(2)}} = {\tau_{2} = {{t_{0}^{(3)} < \ldots < t_{N_{t}}^{(N_{\tau})}} = {\tau_{N_{\tau}} = \Gamma}}}}}}}} \right\}}$

forms a refinement of the partition P_([0, r]) in (6). On the other hand, the space interval [−R, R] can be uniformly discretized by

P _([−R,R]) ={−R=s ₀ <s ₂ <. . . <s _(N) , =R},

wherein s_(j)−s_(j−1)=Δs, j=1, 2, . . . , N_(s) and R is a suitably large number so that the Gaussian distribution can be ignored outside [−R,R]. For the discretization of a D-cell [−R,R]^(D) ⊂ R^(D), an ordered set of the power set of P_([−R,R]), is considered.

P _([−R,R]) ^(D) ={s _(j)}_(j=1) ^((N) ^(s) ⁾ ^(D) ,

wherein s_(j)=(s_(j) ⁽¹⁾, s_(j) ⁽²⁾, . . . , s_(j) ^((D)))^(T), s_(j) ^((d)) ∈ P_([−R,R], j=)1, . . . , (N_(s))^(D), d=1, . . . , D. For the d-th dimension of the space, d=1, . . . , D, the second order partial differential operator can be approximated by using the Euler central difference scheme

$\begin{matrix} {{{\frac{\partial^{2}\overset{\sim}{u}}{\partial s^{2}}\left( {t_{n},s_{j}} \right)} \approx \left\lbrack {{\alpha \left( \frac{U_{j + 1}^{n + 1} - {2U_{j}^{n + 1}} + U_{j - 1}^{n + 1}}{\left( {\Delta \; s} \right)^{2}} \right)} + {\beta \left( \frac{U_{j + 1}^{n} - {2U_{j}^{n}} + U_{j - 1}^{n}}{\left( {\Delta \; s} \right)^{2}} \right)}} \right\rbrack},} & (7) \end{matrix}$

wherein U_(j) ^(n)≡ũ(t_(n),s_(j)) and α+β=1, α,β≧0. Similarly, the partial differential operator can be approximated by

$\begin{matrix} {{\frac{\partial^{2}\overset{\sim}{u}}{\partial s}\left( {t_{n},s_{j}} \right)} \approx {\left\lbrack {{\alpha\left( \frac{U_{j + 1}^{n + 1} - U_{j - 1}^{n + 1}}{2\Delta \; s} \right)} + {\beta \left( \frac{U_{j + 1}^{n} - U_{j - 1}^{n}}{2\Delta \; s} \right)}} \right\rbrack.}} & (8) \end{matrix}$

In other words, the discretized Laplacian operator in (2) can be represented by the matrix

$\begin{matrix} {{T_{d} \equiv \left\lbrack {\left( {\underset{k = 1}{\overset{D - d}{\otimes}}I_{N_{s}}} \right) \otimes T_{d} \otimes \left( {\overset{D}{\underset{k = {D - d + 2}}{\otimes}}{I_{N}}_{s}} \right)} \right\rbrack},} & (9) \end{matrix}$

wherein

denotes the Kronecker product (or tensor product), I_(N) _(s) is the identity matrix of size N_(s) and the matrix

$T_{d} = {{\frac{1}{\left( {\Delta \; s} \right)^{2}}\begin{bmatrix} {- 2} & 1 & \; & \; & \; \\ 1 & {- 2} & 1 & \; & \; \\ \; & 1 & {- 2} & \ddots & \; \\ \; & \; & \ddots & \ddots & 1 \\ \; & \; & \; & 1 & {- 2} \end{bmatrix}}.}$

Similarly, the discretized partial differential operator can be represented by the matrix

$\begin{matrix} {{K_{d} \equiv \left\lbrack {\left( {\underset{k = 1}{\overset{D - d}{\otimes}}{I_{N}}_{s}} \right) \otimes K_{d} \otimes \left( {\overset{D}{\underset{k = {D - d + 2}}{\otimes}}{I_{N}}_{s}} \right)} \right\rbrack},} & (10) \end{matrix}$

wherein the matrix

$K_{d} = {{\frac{1}{2\Delta \; s}\begin{bmatrix} 0 & 1 & \; & \; & \; \\ {- 1} & 0 & 1 & \; & \; \\ \; & {- 1} & 0 & \ddots & \; \\ \; & \; & \ddots & \ddots & 1 \\ \; & \; & \; & {- 1} & 0 \end{bmatrix}}.}$

For each time period t_(n) ^((k)) ∈[τ_(k−1), τ_(k)], the partial differential of time

$\frac{\partial{\overset{\sim}{u}}_{k}}{\partial t}\left( {t_{n},s} \right)$

in (2) can be discretized by

$\begin{matrix} {{{\frac{\partial{\overset{\sim}{u}}_{k}}{\partial t}\left( {t_{n},s} \right)} \approx \frac{U^{{(k)},{n + 1}} - U^{{(k)},n}}{\Delta \; t}},} & (11) \end{matrix}$

wherein U^((k),n)≡(ũ_(k)(t_(n),s₁),ũ_(k)(t_(n), s₂), . . . , ũ_(k)(t_(n), s_((N) _(s) ₎ _(D) ))^(T). Hence the numerical scheme can be written in the form

$\begin{matrix} {{\frac{U^{{(k)},n} - U^{{(k)},{n - 1}}}{\Delta \; t} = {{\alpha \; {AU}^{{(k)},n}} + {\beta \; {AU}^{{(k)},{n - 1}}}}},} & (12) \end{matrix}$

wherein α+β=1, α, β≧0 and the matrix

$\begin{matrix} {{A = {{\frac{1}{2}{\sum\limits_{d = 1}^{D}\; T_{d}}} + {\sum\limits_{d = 1}^{D}\; {P_{d}K_{d}}} + Q}},} & (13) \end{matrix}$

P_(d)≡diag{p_(d)(s₁)}_(j=1) ^((N) ^(s) ⁾ ^(D) , Q≡diag {q(s₁)}_(j=1) ^((N) ^(s) ⁾ ^(D) are diagonal matrices. For each t_(n) ^((k) ∈[τ) _(k−1), τ_(k)], j=1, . . . , N_(t), the linear system is solved

[l−α(Δt) A]U ^((k)n) =[I+β(Δt)A]U ^((k)n−1),   (14)

for k=1, . . . , N_(r) with the initial vector U^((k), 0)≡(U₁ ^((k),0), U₂ ^((k),0), . . . , U_((N) _(s) ₎ _(D) ^((k),0))^(T), in which

$\begin{matrix} {{U_{j}^{{(k)},0} = {\exp \left\{ {\sum\limits_{d = 1}^{M}\; \left\lbrack {{y\left( \tau_{k + 1} \right)} - {{y\left( \tau_{k} \right)}{h_{d}\left( s_{j} \right)}}} \right\rbrack} \right\} U^{{({k - 1})},N,}}},{j = 1},\ldots \mspace{14mu},{\left( N_{s} \right)^{D}.}} & (15) \end{matrix}$

Each vector U^((k),n) in (14) should be normalized such that Σ_(j=1) ^((N) ^(d) ⁾ ^(b) U_(j) ^((k),n)=1. Then the vector U^((k),n) represents the probability distribution of the state at time t_(n) ^((k)). Finally, the expectation is computed

$\begin{matrix} {{\hat{X}\left( t_{n}^{(k)} \right)} = {\sum\limits_{j = 1}^{{(N_{d})}^{D}}\; {s_{j}U_{j}^{{(k)},n}}}} & (16) \end{matrix}$

as the estimation for the real state X(t_(n) ^((k))). In particular, the parameter α=1 and β=0 is chosen since the implicit scheme is stable in most of case while the explicit scheme (α=1, β=1) is usually unstable. The algorithm of “IEM for Nonlinear Filter” for solving the nonlinear filtering problem is stated as the following Algorithm 2.

Input: a terminal time Γ, the space interval [−R, R]^(D), the time step size Δt, the space step size Δs, the vector-valued functions f = (f₁, . . . , f_(D))^(T), h = (h₁, . . . , h_(M))^(T), the observations {Y_(k)}_(k=1) ^(N) ^(t) and the initial state σ₀ Output: the approximation of the state {circumflex over (X)}(t)  1. ${{{Set}\mspace{14mu} {the}\mspace{14mu} {values}\mspace{14mu} N_{\tau}} = {\frac{\Gamma}{\Delta\tau} + 1}},{N_{t} = {{\frac{\Delta\tau}{\Delta t} + {1\mspace{14mu} {and}\mspace{14mu} N_{s}}} = {\frac{2R}{\Delta s} + 1.}}}$  2. for j = 1 to (N_(s))^(D) do  3.   U_(j) ← σ₀(s_(j))exp{Σ_(d=1) ^(M)y_(d)(τ₀)h_(d)(s_(j))}.  4. end for  5. for k = 1 to N_(t) do  6.   for j = 1 to (N_(s))^(D) do  7.   U_(j) ← σ₀(s_(j))exp{Σ_(d=1) ^(M)[y_(d)(τ_(k+1)) − y_(d)(τ_(k))]h_(d)(s_(j))}U_(j) as in (15)  8. end for  9. for n = 1 to N_(t) do 10.   U_(tmp) = [I + β(Δt)A]U. 11.   Solve the linear system [I = β(Δt)A]U = U_(tmp) as in (14) 12.    $\left. {{Normalize}\mspace{14mu} {the}\mspace{14mu} {solution}\mspace{14mu} U}\leftarrow{\frac{U}{\Sigma_{j = 1}^{{(N_{s})}^{D}}U_{j}}.} \right.$ 13.   Set the approximation of state {circumflex over (X)}(t_(n) ^((k))) = Σ_(j=1) ^((Ns)) ^(D) s_(j)U_(j) as in (16) 14.  end for 15. end for The FFTs for the discretized Laplacian matrix Δ is well-known. In the equations (12) and (13), the quasi-implicit scheme is considered as

$\begin{matrix} {{{\left( {1 - {\frac{\Delta \; t}{2}\Delta}} \right)U^{{(k)},n}} = {\left\lbrack {I + {\Delta \; {t\left( {{\sum\limits_{d = 1}^{D}{P_{d}K_{d}}} + Q} \right)}}} \right\rbrack U^{{(k)},{n - 1}}}},} & (17) \end{matrix}$

wherein Δ≡Σ_(d=1) ^(D) T_(d), then the linear system (17) can be efficiently solve by FFTs. For the case of D=1 (one-dimensional case), the Laplacian matrix in (17) satisfies Δ₁≡T₁. By the Fourier sine transformation, the spectral decomposition is:

${T_{1} = {\frac{1}{\left( {\Delta \; s} \right)^{2}}{WSW}^{*}}},$

wherein W≡└W_(ij)┘ with W_(ij)=sin (ijπΔs), and

$s \equiv {{diag}{\left\{ {{- 4}{\sin^{2}\left( \frac{i\; {\pi\Delta}\; s}{2} \right)}} \right\}_{i = 1}^{N_{s}}.}}$

Then the linear system of Δ₁ can be solved by calling the MATLAB function “FFT”. the algorithm of “QIEM for Nonlinear Filter with FFTs” for solving the D-dimensional nonlinear filtering problem by applying the FFT is stated as the following Algorithm 3.

Input: a terminal time Γ, the space interval [−R, R]^(D), the time step size Δt, the space step size Δs, the vector-valued functions f = (f₁, . . . , f_(D))^(T), h = h₁, . . . , h_(M))^(T), the observations {Y_(k)}_(k=1) ^(N) ^(s) and the initial state σ₀ Output: the approximation of the state {circumflex over (X)}(t)  1. ${{{Set}\mspace{14mu} {the}\mspace{14mu} {value}\mspace{14mu} N_{\tau}} = {\frac{\Gamma}{\Delta\tau} + 1}},{N_{t} = {{\frac{\Delta\tau}{\Delta t} + {1\mspace{14mu} {and}\mspace{14mu} N_{s}}} = {\frac{2R}{\Delta s} + 1.}}}$  2. for j = 1 to (N_(s))^(D) do  3.   U_(j) ← σ₀(s_(j))exp{Σ_(d=1) ^(M)y_(d)(τ₀)h_(d)(s_(j))}.  4. end for  5. for k = 1 to N_(τ) do  6.   for j = 1 to (N_(s))^(D) do  7.    U_(j) ← σ₀(s_(j))exp{Σ_(d=1) ^(M)[y_(d)(τ_(k+1)) − y_(d)(τ_(k))h_(d)(s_(j))}U_(j) as in (15)  8.  end for  9.  for n = 1 to N_(t) do 10.   U_(tmp) ← └I + Δt(Σ_(d=1) ^(D)P_(d)K_(d) + Q)┘U. 11.   Call FFTs: U ← (

_(d=1) ^(D)W*)Utmp. 12.   for j = 1 to N_(s) do 13.     $\left. U_{j}\leftarrow{U_{j}/{\left\lbrack {1 + {\frac{2{\Delta t}}{({\Delta s})^{2}}{\sin^{2}\left( \frac{j\; {\pi\Delta s}}{2} \right)}}} \right\rbrack.}} \right.$ 14.   end for 15.   Call IFFTs: U ← (

_(d=1) ^(D)W*)U 16.    $\left. {{Normalize}\mspace{14mu} {the}\mspace{14mu} {solution}\mspace{14mu} U}\leftarrow{\frac{U}{\Sigma_{j = 1}^{{(N_{s})}^{D}}U_{j}}.} \right.$ 17.   Set the approximation of state {circumflex over (X)}(t_(n) ^((k)) = Σ_(j=1) ^((N) ^(s) ⁾ ^(D) s_(j)U_(j) as in (16) 18.  end for 19. end for

The Laplacian matrix in (17) is a second-order approximation of the Laplacian operator. Hereafter, a fourth-order accurate scheme for Laplacian operator which reduces the size of the discretization matrix considerably is considered, but preserves the same accuracy as the second-order approximation. Since there is no general form of the higher-order scheme for Laplacian operator, for convenience in practice, the Kolmogorov equations (2) in two-dimensional case is considered.

The 9-point scheme for the discretized Laplacian operator {tilde over (Δ)}₂ is defined by:

$\begin{matrix} {{{\overset{\sim}{\Delta}}_{2}U_{i,j}} = {\frac{1}{6\left( {\Delta \; s} \right)^{2}}\left\lbrack {{4U_{{i - 1},j}} + {4U_{{i + 1},j}} + {4U_{i,{j - 1}}} + {4U_{i,{j + 1}}} + U_{{i - 1},{j - 1}} + U_{{i + 1},{j - 1}} + U_{{i - 1},{j + 1}} + U_{{i + 1},{j + 1}} - {20U_{i,j}}} \right\rbrack}} & (19) \end{matrix}$

which is a fourth-order approximation of Laplacian operator. The matrix form of (19) is represented as

$\begin{matrix} {{{{\overset{\sim}{\Delta}}_{2} = {\frac{1}{\left( {\Delta \; s} \right)^{2}}\begin{bmatrix} \Sigma & \Phi & \; & \; \\ \Phi & \Sigma & \ddots & \; \\ \; & \ddots & \ddots & \Phi \\ \; & \; & \Phi & \Sigma \end{bmatrix}}},{wherein}}{\Sigma = {\frac{1}{3}\begin{bmatrix} {- 10} & 2 & \; & \; \\ 2 & {- 10} & \ddots & \; \\ \; & \ddots & \ddots & 2 \\ \; & \; & 2 & {- 10} \end{bmatrix}}},{\Phi = {{\frac{1}{6}\begin{bmatrix} 4 & 1 & \; & \; \\ 1 & 4 & \ddots & \; \\ \; & \ddots & \ddots & 1 \\ \; & \; & 1 & 4 \end{bmatrix}}.}}} & (20) \end{matrix}$

Since the matrix {tilde over (Δ)}₂ is also a Toeplitz matrix, the fast Fourier transformation is derived for solving the linear system {tilde over (Δ)}₂U=b. Note that

$\begin{matrix} {{\overset{\sim}{\Delta}}_{2} = \left( {{\frac{1}{6}\left\lbrack {\left( {T_{1} + {6I}} \right) \otimes \left( {T_{1} + {6I}} \right)} \right\rbrack} - {6I}} \right)} \\ {= {\frac{1}{6\left( {\Delta \; s} \right)^{2}}\left( {\left\lbrack {\left( {{WSW}^{*} + {6I}} \right) \otimes \left( {{WSW}^{*} + {6I}} \right)} \right\rbrack - {36I}} \right)}} \\ {= {\frac{1}{6\left( {\Delta \; s} \right)^{2}}\left( {{\left( {W \otimes W} \right)\left( {\left( {S + {6I}} \right){W^{*} \otimes \left( {S + {6I}} \right)}W^{*}} \right)} - {36I}} \right)}} \\ {= {\frac{1}{6\left( {\Delta \; s} \right)^{2}}\left( {\left( {W \otimes W} \right)\left( {\left( {\left( {S + {6I}} \right) \otimes \left( {S + {6I}} \right)} \right) - {36\left( {I \otimes I} \right)}} \right)\left( {W^{*} \otimes W^{*}} \right)} \right)}} \end{matrix}$

wherein

$T_{1} = {\frac{1}{\left( {\Delta \; s} \right)^{2}}{WSW}^{*}}$

as given in (18). Based on the QIEM Algorithm 3, the algorithm “Fourth-order QIEM for 2-D Nonlinear Filter with FFTs” for solving the two-dimensional nonlinear filtering problem by applying the fourth-order QIEM with FFTs is stated in the following Algorithm 4.

Input: a terminal time Γ, the space interval [−R, R], the time step size Δt, the space step size Δs, the functions f = (f₁, f₂)^(T), h = (h₁, h₂)^(T), the observations {Y(τ_(k)) = (y₁(τ_(k)), y₂(τ_(k)))^(T)}_(k=1) ^(N) ^(s) and the initial state σ₀ Output: the approximation of the state {circumflex over (X)}(t)  1. ${{{Set}\mspace{14mu} {the}\mspace{14mu} {values}\mspace{14mu} N_{\tau}} = {\frac{\Gamma}{\Delta\tau} + 1}},{N_{t} = {{\frac{\Delta\tau}{\Delta t} + {1\mspace{14mu} {and}\mspace{14mu} N_{s}}} = {\frac{2R}{\Delta s} + 1.}}}$  2. for j = 1 to (N_(s))² do  3.  U_(j) ← σ₀(s_(j))exp{y₁(τ₀)h₁(s_(j)) + y₂(τ₀)h₂(s_(j))}.  4. end for  5. for k = 1 to N_(τ) do  6.  for j = 1 to N_(s) do  7.   U_(j) ← exp{[y₁(τ_(k+1)) − y₁(τ_(k))]h₁(s_(j)) + [y₂τ_(k+1)) − y₂(τ_(k))]h₂(s_(j))}U_(j)  8.  end for  9.  for n = 1 to N_(t) do 10.   U_(imp) ← [I + Δt(P₁D₁ + P₂D₂ + Q)]U. 11.   % Here the matrixes W and S are defined in (18). 12.   Call FFTs: U ← (W* 

 W*)U_(imp). 13.   for j = 1 to N_(s) do 14. $\left. U_{j}\leftarrow{\left\lbrack {\left( {I \otimes I} \right) - {\frac{\Delta t}{12({\Delta s})^{2}}\left( {{\left( {S + {6I}} \right) \otimes \left( {S + {6I}} \right)} - {36\left( {I \otimes I} \right)}} \right)}} \right\rbrack^{- 1}{U_{j}.}} \right.$ 15.   end for 16.   Call FFTs: U ← (W 

 W)U. 17.    $\left. {{Normalize}\mspace{14mu} {the}\mspace{14mu} {solution}\mspace{14mu} U}\leftarrow{\frac{U}{\Sigma_{j = 1}^{N_{s}}U_{j}}.} \right.$ 18.   Set the approximation of state {circumflex over (X)}(t_(n) ^((k))) = Σ_(j=1) ^((N) ^(s) ⁾ ² s_(j)U_(j) 19.  end for 20. end for The computation of nonlinear filtering problem is a real time problem. Saving the computational cost becomes an essential issue. In order to solve the nonlinear filtering problem in a more efficient way, the superposition technique is adopted, and the Dirac delta functions is

{δ_(C) _(i) (S)=e ^(−η|s−C) ^(k) ^(|) ¹ |C _(k)=(c _(k) _(l) , . . . , c _(k) _(D) )^(τ) ∈[−R, R] ^(D)}_(k=i) ^(N) ^(s)

with S=(s_(i), . . . , s_(D))^(τ) ∈R^(D) and ∥S−C_(k)∥²=Σ_(j=1) ^(D)(s_(j)−c_(k) _(j) )², as various initials to compute the approximate states for the nonlinear filtering problem, separately. Then all the fundamental solutions {ν_(k)}_(k=1) ^(N) ^(s) corresponding to the Dirac delta functions {δ_(c) _(k) }_(k=1) ^(N) ^(D) are stored.

In practice, for any given initial probability density function u₀, a set of coefficients {α_(k)}_(k=1) ^(N) ^(s) of the linear combination of Dirac delta functions is calculated to satisfy

$u_{0} \approx {\sum\limits_{k = 1}^{N_{s}}{\alpha_{k}{\delta_{C_{k}}.}}}$

Then the approximate probability density function of the state v can be directly obtained by computing the linear combination of the fundamental solutions

$v \approx {\sum\limits_{k = 1}^{N_{s}}{\alpha_{k}{v_{k}.}}}$

This method significantly saves a large amount of computational cost.

The linear operator defined in (14) is a nonnegative operator, the solution of each time step represents a probability distribution of the space. In order to guarantee the property that each solution is nonnegative, it is found that the sufficient condition such that the matrix (I−ΔtA)⁻¹ is a nonnegative operator. First, the definition of an M-matrix and its equivalence condition are as follows.

Definition: A real matrix B=└B_(ij)┘ is called an M-matrix if B_(ij)≦0, i≠j and B⁻¹ exists with B⁻¹≧0. Lemma : (Equivalence Condition of M-matrix). Let B be a real matrix with B_(ij)≦0 for i≠j. Then B is an M-matrix if and only if there is a positive vector ν>0 such that B_(ν)>0. The following thereon shows that the vector U in each iteration of step 11 in Algorithm 2 preserves nonnegativity of the probability density function. Theorem 1: (Sufficient Condition of Nonnegative Operator). Given real-valued functions p_(d), d=1, . . . , D, q, a time step Δt and a space step Δs. Let B≡I−ΔtA, wherein A (13). If for each s ∈[−R, R]^(D),

$\begin{matrix} {{{{p_{d}(S)}} < \frac{1}{\Delta \; s}},{{{q(S)}} < \frac{1}{\Delta \; t}},} & (21) \end{matrix}$

for d=1, . . . , D, the B is an M-matrix. That is, B⁻¹≧0 is a nonnegative operator. Proof. First to check B_(ij)≦0 for i≠j. From the structure of the matrices in (9) and (10) that for i≠j either B_(ij)=0 or

$\begin{matrix} \begin{matrix} {B_{ij} = {{- \Delta}\; {t\left( {\frac{1}{2\left( {\Delta \; s} \right)^{2}} + \frac{p_{d}(S)}{2\Delta \; s}} \right)}}} \\ {= {{- \frac{\Delta \; t}{2\left( {\Delta \; s} \right)^{2}}}\left( {1 + {\Delta \; {{sp}_{d}(S)}}} \right)}} \\ {\leq {{- \frac{\Delta \; t}{2\left( {\Delta \; s} \right)^{2}}}\left( {1 - {\Delta \; s{{p_{d}(S)}}}} \right)}} \\ {< 0} \end{matrix} & (22) \end{matrix}$

Then inequality (22) follows from the first equation of (21). Next, B1>0 is checked, wherein 1≡(1, 1, . . . , 1)^(τ)>0. Note that B1 is a vector whose entry is the row sum of B. Hence

$\begin{matrix} \begin{matrix} {{B\; 1} = {1 - {\Delta \; {t\left( {\frac{- k}{2\left( {\Delta \; s} \right)^{2}} + \frac{{\Sigma_{d = 1}^{k}\left( {- 1} \right)}^{m_{d}}{p_{d}(S)}}{2\Delta \; s} + {q(S)}} \right)}}}} \\ {\geq {1 - {\Delta \; {t\left( {\frac{- k}{2\left( {\Delta \; s} \right)^{2}} + \frac{{kmax}_{d}{{p_{d}(S)}}}{2\Delta \; s} + {q(S)}} \right)}}}} \\ {= {\left( {1 - {\Delta \; {{tq}(S)}}} \right) + {\frac{k\; \Delta \; t}{2\left( {\Delta \; s} \right)^{2}}\left( {1 - {\Delta \; s\; {\max\limits_{d}{{p_{d}(S)}}}}} \right)}}} \\ {\geq {\left( {1 - {\Delta \; t{{q(S)}}}} \right) + {\frac{k\; \Delta \; t}{2\left( {\Delta \; s} \right)^{2}}\left( {1 - {\Delta \; s\; {\max\limits_{d}{{p_{d}(S)}}}}} \right)}}} \\ {{> 0},} \end{matrix} & (23) \end{matrix}$

for some k ∈{1, . . . , D}, m_(d) ∈{0, 1}. The inequality (23) follows from (21). By Lemma 1, B is an M-matrix. That is, B⁻¹≧0. Consequently, by Theorem 1, the vector U=[I−ΔtA]⁻¹U_(tmp) in Step 11 of Algorithm 2 is nonnegative.

The convergence of the IEM and the QIEM is proved by checking the consistency and the stability of the schemes.

Theorem 2: (Consistency of IEM.QIEM). The local truncation errors of IEM (12) and QIEM (17) are O(Δt+(Δs)²). That is, IEM and QIEM are consistent. Proof. The first-order Taylor expansion of u at the point (t+Δt,s) implies

$\begin{matrix} {{\frac{\partial u}{\partial t}\left( {t,s} \right)} = {\frac{{u\left( {{t + {\Delta \; t}},s} \right)} - {u\left( {t,s} \right)}}{\Delta \; t} + {{O\left( {\Delta \; t} \right)}.}}} & (24) \end{matrix}$

The third-order Taylor expansions of u at the points (t,s+Δs) and (t,s−Δs), respectively, lead to

$\begin{matrix} {{u\left( {t,{s + {\Delta \; s}}} \right)} = {{u\left( {t,s} \right)} + {\Delta \; s\frac{\partial u}{\partial s}\left( {t,s} \right)} + {\frac{\left( {\Delta \; s} \right)^{2}}{2}\frac{\partial^{2}u}{\partial s^{2}}\left( {t,s} \right)} + {\frac{\left( {\Delta \; s} \right)^{3}}{6}\frac{\partial^{3}u}{\partial s^{3}}\left( {t,s} \right)} + {O\left( \left( {\Delta \; s} \right)^{4} \right)}}} & (25) \\ {\mspace{79mu} {and}} & \; \\ {{u\left( {t,{s - {\Delta \; s}}} \right)} = {{u\left( {t,s} \right)} + {\Delta \; s\frac{\partial u}{\partial s}\left( {t,s} \right)} + {\frac{\left( {\Delta \; s} \right)^{2}}{2}\frac{\partial^{2}u}{\partial s^{2}}\left( {t,s} \right)} - {\frac{\left( {\Delta \; s} \right)^{3}}{6}\frac{\partial^{3}u}{\partial s^{3}}\left( {t,s} \right)} + {{O\left( \left( {\Delta \; s} \right)^{4} \right)}.}}} & (26) \end{matrix}$

By adding the equations (25) and (26), obtaining

$\begin{matrix} {{\frac{\partial^{2}u}{\partial s^{2}}\left( {t,s} \right)} = {\frac{{u\left( {t,{s + {\Delta \; s}}} \right)} - {2{u\left( {t,s} \right)}} + {u\left( {t,{s - {\Delta \; s}}} \right)}}{\left( {\Delta \; s} \right)^{2}} + {{O\left( \left( {\Delta \; s} \right)^{2} \right)}.}}} & (27) \end{matrix}$

Similarly, by subtracting the equation (26) rom (25), obtaining

$\begin{matrix} {{\frac{\partial u}{\partial s}\left( {t,s} \right)} = {\frac{{u\left( {t,{s + {\Delta \; s}}} \right)} - {u\left( {t,{s - {\Delta \; s}}} \right)}}{2\Delta} + {{O\left( \left( {\Delta \; s} \right)^{2} \right)}.}}} & (28) \end{matrix}$

Hence, according to the equations (7), (8) and (11), respectively, the equations (27), (28) and (24) show the local truncation error of (12) is O(Δt+(Δs)²). Theorem 3: (Sufficient Condition for Stability of IEM). The IEM (12) is stable if the function f in (1) satisfies

∇·ƒ≧0.   (29)

Proof. IEM (12) is stable by applying von Neumann stability analysis. Let U_(j) ^(n)=ξ(k)^(n)e^(tkj(Δs)), wherein t≡√{square root over (−1)} and ξ(k) is known as the amplification factor. Substituting U_(j) ^(n) into the scheme (12), obtaining

$\frac{{\xi (k)} - 1}{\Delta \; t} = {{\frac{\xi (k)}{2\left( {\Delta \; s} \right)^{2}}\left( {e^{{ikj}{({\Delta \; s})}} - 2 + e^{- {{ikj}{({\Delta \; s})}}}} \right)} + {\frac{\xi (k)}{2\Delta \; s}\left( {e^{{ikj}{({\Delta \; s})}} - e^{{ikj}{({\Delta \; s})}}} \right){p(x)}} + {{\xi (k)}{q(s)}}}$

That is,

$\begin{matrix} {{\xi (k)} = \left\lbrack {1 - {\Delta \; {t\left( {\frac{e^{{ikj}{({\Delta \; s})}} - 2 + e^{{ikj}{({\Delta \; s})}}}{2\left( {\Delta \; s} \right)^{2}} + {\frac{e^{{ikj}{({\Delta \; s})}} - e^{{ikj}{({\Delta \; s})}}}{2\Delta \; s}{p(s)}} + {q(s)}} \right)}}} \right\rbrack^{- 1}} \\ {= \left\lbrack {1 - {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {{\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - 1 + {i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s} + {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)}} \right\rbrack^{- 1}} \\ {= \begin{bmatrix} {1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)} -} \\ {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s} \end{bmatrix}^{- 1}} \\ {= {\frac{1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)} + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s}}{\begin{matrix} {\left( {1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)}} \right)^{2} +} \\ \left( {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s} \right)^{2} \end{matrix}}.}} \end{matrix}\quad$

If ∇·ƒ≧0, then

${q(s)} \equiv {- \left\lbrack {{\nabla{\cdot f}} + {\frac{1}{2}{\sum\limits_{j = 1}^{M}{h_{j}^{2}(S)}}}} \right\rbrack} \leq 0.$

It follows that

$\begin{matrix} {{{\xi (k)}}^{2} = \frac{\begin{matrix} {\left( {1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)}} \right)^{2} +} \\ \left( {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s} \right)^{2} \end{matrix}}{\begin{pmatrix} {\left( {1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)}} \right)^{2} +} \\ \left( {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}\Delta \; s} \right)^{2} \end{pmatrix}^{2}}} \\ {= \begin{bmatrix} {\left( {1 + {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {1 - {\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - {{q(s)}\left( {\Delta \; s} \right)^{2}}} \right)}} \right)^{2} +} \\ \left( {\frac{\Delta \; t}{\Delta \; s}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}} \right)^{2} \end{bmatrix}^{- 1}} \\ {\leq \left\lbrack \left( {1 - {\Delta \; {{tq}(s)}}} \right)^{2} \right\rbrack^{- 1} \leq 1} \end{matrix}\quad$

This implies that IEM (12) is stable under the assumption that ∇∫ƒ≧0. Theorem 4: (Sufficient Condition for Stability of OIEM). The QIEM (17) is stable if both the step size of time Δt and the step size of space Δs are sufficient small. More precisely, ≢t and Δs satisfy

(Δs)²(2q(s)+q(s)² Δt)+Δtp(s)²≦2.   (30)

Proof. As in the proof of Theorem 3, U_(j) ^(N)=ξ(k)^(n) e^(tkj(Δt)) is substituted into the scheme (17) to obtain

$\frac{{\xi (k)} - 1}{\Delta \; t} = {{\frac{\xi (k)}{2\left( {\Delta \; s} \right)^{2}}\left( {e^{{ikj}{({\Delta \; s})}} - 2 + e^{{ikj}{({\Delta \; s})}}} \right)} + {\frac{1}{{2\Delta \; s}\;}\left( {e^{{ikj}{({\Delta \; s})}} - e^{{ikj}{({\Delta \; s})}}} \right){p(s)}} + {{q(s)}.}}$

That is,

${\xi (k)} = {\frac{1 + {\Delta \; {{tq}(s)}} + {\frac{\Delta \; t}{\Delta \; s}i\; {\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}}}{1 - {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {{\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - 1} \right.}}.}$

If Δs and Δt satisfy (Δs)²(2q(s)²Δt)+Δtp(s)²≦2, then

${{\left( {\Delta \; s} \right)^{2}\left( {{2{q(s)}} + {{q(s)}^{2}\Delta \; t}} \right)} + {\Delta \; {{tp}(s)}^{2}}} \leq {2 + {\frac{{\Delta \; t}\;}{\left( {\Delta \; s} \right)^{2}}.}}$

Multiplying both side by

$\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}},$

having

${{2{q(s)}\Delta \; t} + {{q(s)}^{2}\left( {\Delta \; t} \right)^{2}} + {\frac{\left( {\Delta \; t} \right)^{2}}{\left( {\Delta \; s} \right)^{2}}{p(s)}^{2}}} \leq {\frac{2\Delta \; t}{\left( {\Delta \; s} \right)^{2}} + {\frac{\left( {\Delta \; t} \right)^{2}}{\left( {\Delta \; s} \right)^{4}}.}}$

Adding both side by 1, obtaining

${\left( {1 + {\Delta \; {{tq}(s)}}} \right)^{2} + \left( {\frac{\Delta \; t}{\Delta \; s}{p(s)}} \right)^{2}} \leq {\left( {1 + \frac{\Delta \; t}{\left( {\Delta \; s} \right)}} \right)^{2}.}$

It follows that

${{\xi (k)}}^{2} = {\frac{\left( {1 + {\Delta \; {{tq}(s)}}} \right)^{2}\left( {\frac{\Delta \; t}{\Delta \; s}{\sin \left( {{kj}\left( {\Delta \; s} \right)} \right)}{p(s)}} \right)^{2}}{\left( {1 - {\frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}\left( {{\cos \left( {{kj}\left( {\Delta \; s} \right)} \right)} - 1} \right)}} \right)^{2}} \leq \frac{\left( {1 + {\Delta \; {{tq}(s)}}} \right)^{2} + \left( {\frac{\Delta \; t}{\Delta \; s}{p(s)}} \right)^{2}}{\left( {1 + \frac{\Delta \; t}{\left( {\Delta \; s} \right)^{2}}} \right)^{2}} \leq 1}$

This implies QIEM (17) is stable under the given assumption. Theorem 5: (Sufficient Conditions for Convergence of IEM/QIEM). The IEM and QIEM converge if the conditions of (29) and (30) hold, respectively. Proof. From the consistency of IEM/QIEM in Theorem 2 as well as the stabilities of IEM and QIEM in Theorem 3 and Theorem 4, respectively, the convergence if IEM/QIEM follows by the Lax-Richtmyer equivalence theorem immediately.

In summary, the method for solving high-dimensional nonlinear filtering problems of the preset invention has the following advantages compared with algorithms and schemes available now.

1. The method of the present invention solves D-dimensional nonlinear filtering problems and gets approximate numerical solutions based on Yau-Yau filtering theory. The Quasi-Implicit Euler Method (QIEM) is applied to solve the Kolmogorov equations and estimate approximate numerical solutions of the signal-observation model. Moreover, QIEM is feasible for acceleration by fast Fourier transformation (FFT). Thus stability of the numerical solutions is ensured and a large amount of computational cost is saved.

2. The method of the present invention guarantees nonnegativity of the numerical solutions of Kolmogorov equations in each iteration and preserves probability density functions in the iterative process. The numerical results show that the method is efficient and promising.

Additional advantages and modifications will readily occur to those skilled in the art. Therefore, the invention in its broader aspects is not limited to the specific details, and representative devices shown and described herein. Accordingly, various modifications may be made without departing from the spirit or scope of the general inventive concept as defined by the appended claims and their equivalents. 

What is claimed its:
 1. A method for solving high-dimensional nonlinear filtering problems comprising the steps of: solving an equation by a fast computational module to obtain approximate numerical solutions of a signal-observation model; and accelerating a process of solving the equation by the fast computational module for improving computational stability.
 2. The method claimed in claim 1, wherein a Quasi-Implicit Euler Method(QIEM) is applied to solve Kolmogorov equations and obtain the approximate numerical solutions of the signal-observation model.
 3. The method as claimed in claim 2, wherein the Quasi-Implicit Euler Method (QIEM) is iteratively formulated by: ${{\left\lbrack {I_{N^{D}} - {\frac{\Delta \; t}{2\left( {\Delta \; s} \right)^{2}}L_{N}^{(D)}}} \right\rbrack u_{n + 1}} = {\left\lbrack {I_{N^{D}} + {\Delta \; {t\left( {{\frac{1}{2\Delta \; s}K_{\; N}^{(D)}} + Q_{N}^{(D)}} \right)}}} \right\rbrack u_{N}}};$ wherein ${L_{N}^{(D)} = {\sum\limits_{d = 1}^{D}\left\lbrack {I_{N^{D - d}} \otimes L_{N} \otimes I_{N^{d - 1}}} \right\rbrack}},{{K_{N}^{(D)} = {\sum\limits_{d = 1}^{D}{P_{d}\left\lbrack {I_{N^{D - d}} \otimes K_{N} \otimes I_{N^{d - 1}}} \right\rbrack}}};}$ wherein I_(N) is the identity matrix of size N and P_(d)=diag{p_(d)(s₁)}_(j=1) ^(N) ^(D) , Q=diag{q(s₁)}_(j=1) ^(N) ^(D) ; wherein the matrices L_(N) and K_(N) are defined by ${L_{N} = \begin{bmatrix} {- 2} & 1 & \; & \; \\ 1 & {- 2} & \ldots & \; \\ \; & \ldots & \ldots & 1 \\ \; & \; & 1 & {- 2} \end{bmatrix}},{K_{N} = {\begin{bmatrix} 0 & 1 & \; & \; \\ {- 1} & 0 & \ldots & \; \\ \; & \ldots & \ldots & 1 \\ \; & \; & {- 1} & 0 \end{bmatrix}.}}$ 